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Abstract 

Exact solution to the hierarchy of nonlinear lattice Boltzmann (LB) kinetic equations in the 
stationary planar Couette flow is found at non-vanishing Knudsen numbers. A new method of 
solving LB kinetic equations which combines the method of moments with boundary conditions 
for populations enables to derive closed- form solutions for all higher-order moments. Convergence 
of results suggests that the LB hierarchy with larger velocity sets is the novel way to approximate 
kinetic theory. 
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y important 
l|. Recently, 



Emerging field of fluid dynamics at a micrometer scale becomes increasim 
due to fundamental engineering issues of micro-electromechanical systems 
much attention was attracted by the use the lattice Boltzmann (LB) models for simulation 

8|. By now, it is understood that 
ry based on discrete velocity sets with 
or rational-number approximations 



of microflows by a number of groups 
lattice Boltzmann models form a well-defined hierarc 
velocities defined as zeroes of Hermite polynomials 



thereof 



101 ] . The LB hierarchy constitutes a novel approximation of the Boltzmann equation 



and has to be considered as an alternative to more standard approaches such as higher-order 
hydrodynamics (Burnett or super-Burnett) or Grad's moment systems (for a review see, e. 
g. III). One salient feature of the LB hierarchy, which is crucial to the present study 
and eventually to any realistic application, is that it is naturally equipped with relevant 
boundary conditions derived from Maxwell-Boltzmann's theory 2|. 

Agreement between the LB simulations and kinetic theory j^j, hydrodynamics with slip 
boundary conditions and molecular dynamics [8[| was reported. However, most of these 
numerical works rely on simulation with finite accuracy while the crucial question whether 
or not the kinetic equations underpinning the LBM method are valid physical models of 
microflow remains unanswered. Therefore, it is not surprising to read comments claiming, 
for example, that the slip flow in the LBM is due to discretization errors rather than a 
physical effect, ([3] and references therein). Thus, validity of LBM cannot be addressed 
unless a comparison to representative exact solutions is performed. It is needless to say that 
exact solutions to nonlinear kinetic equations in realistic geometries are very rare. 

In this Letter, we show that the LB hierarchy of kinetic models admits a much more 
accurate treatment. In particular, we find closed- form analytical solutions to nonlinear 
kinetic equations of the LB hierarchy in the stationary planar Couette flow. Not only the slip 
velocity, but also the shear stress and the normal stress difference are evaluated in a closed 
form. Comparison to the kinetic theory demonstrates convergence of approximations with 
the increase of the number of velocities. In the nonlinear domain, even the first member of the 
LB hierarchy predicts nontrivial normal stress which is confirmed with a more microscopic 
direct simulation Monte Carlo (DSMC) method jl3|. The accurate results obtained herein 
strongly suggest that the LB hierarchy should be considered as a novel general tool of kinetic 
theory rather than a plain solver for hydrodynamics. 

Kinetic equations studied in this paper are two-dimensional isothermal discrete velocity 
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models with the Bhatnagar-Gross-Krook (BGK) nonlinear collision integral (for a derivation 
of these models from the Boltzmann-BGK equation see, e. g., 

,Q), 

d t fi + c ia d a f i ^--(f i -fP), (1) 

T 

where summation convention is applied, r is the relaxation time, and the equilibrium dis- 
tribution is 

fP{p, jxJv) = w i yp + 3 ~^r- + j0 - c 2 s 5 a(3 ) J . (2) 

Here c s = a/ (k B T )/m is the speed of sound, T is the reference temperature. 

The first member of the LB hierarchy is the so-called D2Q9 model where the discrete 
velocities Cj a , i — 0, ... ,8, and the weights Wi are 

c, = V3c s {0, 1, 0, -1, 0, +1, -1, -1, +1} 

c y = V3c s {0, 0, 1, 0, -1, +1, +1, -1, -1} (3) 
w = (1/36) {16, 4, 4, 4, 4, 1,1, 1,1}. 

The model (CQ) conserves locally the density p = Yli=o & anc ^ ^ ne momentum density j a = 
Y^i=o c iafi but not the energy. In the hydrodynamic limit, it recovers the Navier-Stokes 
equations with the kinematic viscosity v = tc 2 . 

We consider the planar Couette flow, where a fluid is enclosed between two parallel plates 
separated by a distance L. The bottom plate at y — —L/2 moves with the velocity U\ and 
top plate at y = L/2 moves with the velocity U2 (see Fig. [1]). Let us introduce the mean 
free path / = \/3tc s and the Knudsen number Kn = l/L. The solution for the x- velocity of 
the D2Q9 model derived below reads: 




FIG. 1: Couette flow geometry. Discrete velocities of the D2Q9 model at the bottom and the top 
plates are indicated to explain boundary conditions. 
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where 

6 9 = 1 + 2Kn, (5) 

and where AU = U 2 — U\ is the relative velocity of the plates, and U = (U\ + L^)/2 is the 
centerline velocity. Solution (jlj) is Galilean invariant: if a constant velocity C is added to 
both the plates, U[ 2 = C/i^ + C then the velocity (BJ transforms as = m x + C. The slip 
velocity at the plates, u x (±L/2), features the expected behavior from a linear increase with 
Kn at small Kn to a plug flow at Kn — > oo where u x becomes position- independent. Recent 
careful numerical study of the D2Q9 model by Sofonea and Sekerka [jj revealed the same 
result (gD and ©. 

The next member of the LB hierarchy is the model based on the roots of fourth-order 
Hermite polynomial {±a, ±6}, where a = a/3 — a/6 and b = a/3 + a/6. In two dimensions, 
the discrete velocities are all possible tensor products of the two copies of the four-sets 



{±o, ±6}. For this D2Q16 model [<i 



[lo| . the solution for the velocity profile is found to be: 



where 



»*=ih sinh (ih) AU+ ir e (i) AU+u - (6) 



^^(^k) + 2 v / 3sinh(^) i 



'16 



(4Kn + fj) cosh C^^j + 2 (//Kn + v^) sinh f J J 



(8) 



4Kn 

and /i = a + 6 3.076. The difference between (j4| and (jSJ) is that the latter predicts the 
boundary Knudsen layer (first term in (jHJ)) in a qualitative agreement with kinetic theory 
Jjj]. On the quantitative side, our analytical results can be immediately compared with 



the classical study of the linearized Boltzmann-BGK equation by Willis [15| (see Fig. [2, 
where also the results of the DSMC simulation are reported; note that data in Fig. [2] is 
parameterized with the Knudsen number according to a relation, Kn = where 



a = — yj 2k^T ^ s a P arame ter used in Table I of Ref. [15j). While the simplest D2Q9 model 
predicts well a slip- flow jso 
with numerical studies 



ution, it fails in the transient regime (Kn > 0.1), in agreement 



16]. However, already the D2Q1Q model considerably improves 
the situation. The strong pattern of convergence with increasing the number of velocities in 
the LB hierarchy is clearly there. 
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FIG. 2: Comparison of the LB hierarchy with the linearized Boltzmann-BGK model [15I ] and 
DSMC simulation. Left: Slip velocity at the wall as a function of Knudsen number. Right: Slope 
of the velocity profile at the centerline. Plotted is the deviation from the Navier-Stokes prediction, 
Y = AU- 1 (du x /dy)\ y=Q -l. 

We shall now proceed with major steps of derivation for the D2Q9 model. Firstly, we write 
the kinetic equation for nine populations (PQ) in a form of a moment system for nine moments 
which we choose as follows: Three locally conserved fields, p, j x , j y ; three independent 
components of the pressure tensor, P af3 = ^i=o /« Cia c *^' wmcrl we choose as the trace 
P = Pxx + Pyy, N = P xx — P yy (normal stress difference), and P xy ; two components of the 
energy flux, q a = Y^=o fi°ai c h anc ^ a scalar fourth-order moment, ip = R y yyy+R xxxx —2R XX y yj 
where R a ^e = Y^i=o f iCiaCi P Ci i CiS - ^ ne resulting moment system includes three balance 
equations for the locally conserved fields, 



d t p + d x j x + dyjy 
dt3x + \dx{P + N)+d y P xy 

d t Jy + d x P xy + ^dy(P-N) 



0. 
0. 



(9) 
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three equations for the components of the pressure tensor, 

d t Pxy+d x (q y - 3c 2 s j y ) + d y (q x - 3c 2 j x ) = _ p xy ^ 



r p 

d t N+d x (6c 2 Jx - q x ) + d y (q y - 6c 2 j y ) = i - ^ - N), 



d t P+d x q x + d y q y = ±(2pc 2 s + J j- P), (10) 

two equations for the components of the energy flux, 

d t q x + d x [3c 2 s (P + \n) - + Qc 2 d y P xy = - (Ac 2 Jx - q x ), 

d t q y + 6c 2 s d x P xy + d y [3c 2 s (P - \n) - U) = -(4c 2 s3y - q y ), (11) 

Z Z T 

and, finally, the equation for the fourth-order moment ij), 

d t ^ + 3 c 2 s d x (6c 2 3x -q x ) + 3 c 2 s d y (6c 2 sh - q y ) = \a P c a s + cf- - (12) 

Secondly, we find steady state solution to the moment system under two conditions: (i) 
Unidirectional flow: As the plates extend to infinity in the x direction, we can expect that 
the steady state solution will be independent of x, and (ii) Impermeable plates: The normal 
mass flux equals zero at the walls. For a unidirectional stationary flow, balance equations 
([9]) read: d y j y = 0, d y P xy = 0, and d y (P — N) = 0, whereupon, using condition (ii), we 
get j y = 0, P xy = P X y q , P = N + P , where integration constants P™ q and Po will be 
determined below (superscript 'neq' is added to emphasize that only the non-equilibrium 
part is nontrivial for P xy ). Furthermore, equation for the pressure tensor (TTOT) reads 

d y (q x -3c 2 sJx ) = ~P^, (13) 
d y q v = -{^-N) } (14) 

T p 

d y q y = - (2pc 2 s + J -^-N-P ), (15) 

T P 

From (fl4"I) and (|T5l) it follows Po = 2pc 2 . Thus, the stationary density is constant. Equation 
for energy flux ffTT]) simplifies as: 

q x = 4c 2 s j x , (16) 
d y [5c 2 AP-\N)-^} = - l -q y . (17) 
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Substituting (TT6l) into ({13]) . and integrating the resulting differential equation for j x , we 
obtain the result for the nontrivial velocity component u x = j x /p, 

u x {y) = -— 2 {y + rV) 1 (18) 

T C s 

where V is constant of integration, and where we have introduced tt = P X y q / p. Thus, we find 
the solution for the velocity (up to two constants, tt and V) before higher-order moments 
are addressed. We note in passing that it is precisely the relation ffT6]) pertinent to the 
low-symmetry D2Q9 model (the energy flux is proportional to the momentum flux) which 
precludes the development of the boundary Knudsen layer. This constraint is removed in 
the more symmetric D2Q16 model and in any higher-order member of the LB hierarchy. 

The stationary equation for fourth order moment ip ( Tl2|) . together with ( [Ml) , gives ip = 
Ac 2 s & - 3c 2 s N + 4pcf. Finally, from ([171) and (ED, we get 



iv a Arneq _ y 3/ P 



9ylT q = -— 9y NUeq = - fviO- ( 19 ) 

The ordinary differential equations f|T9|) can be integrated explicitly since the velocity u x (y) 
is available from (Tl8l) . Denoting (p(y) = exp (y/\/3rc s ), the result is 

V3c s N^ + q^ = Acp(-y)-^-\y + T(V-V3c s )(l-<p(-y)) , 

V3c s N ne * - qp = B<p{y) + ^ \y + r(V + V3c s ) (1 - cp( y ))} , (20) 

rcf L J 

where A and B are constants of integration. Thus, the solution to the stationary moment 

system depends on the four integration constants, tt, V, A and B. In order to determine 

these, we need to specify boundary conditions at the moving plates. Note that this is 

precisely where the LB hierarchy differs from the method of moments. It is well known 

that for moment methods, such as Grad systems, it is not possible to provide self-consistent 

boundary conditions for the moments. In our case, this is possible because the boundary 

conditions for the LB equations are formulated in terms of populations rather than in terms 

of moments. Upon inverting the linear relation between the moments and the populations, 

and using the solution for the moments derived above, we obtain the stationary populations 

ft = fP + ,A nec \ where the stationary equilibrium part is given by ([2]) with j y = and 

j x = pu x ( TTBl . while the non-equilibrium part has the form, 

/ Pncq „neq /Vneq \ 

rneq _ I ± xy ( 2 _/l 2 "\ _i_ ( 2 _ 2 "\ 2 I (9~\ \ 

Ji — Uli I ^ 4 Ci x Ci y -\- \ c iy c i ^ c s c iy)< \ ix Cs / C W J ' V^^J 



Thirdly and finally, we apply the classical diffuse boundary conditions [14], which were 
adapted to the present model in Ref. At the bottom plate (y = —L/2), diffuse boundary 
condition in the steady state is (see Fig. [Q, 

/2,5,6|, = _ L/2 =/ 2 ° q 5, 6 (p,pf/i,0). (22) 

In other words, in the steady state, the diffuse boundary condition reduces to setting the 
corresponding populations at equilibrium (j2j) with the density p and velocity of the wall. 
Now, in order to find a relation between the two integration constants, V and tt, we no- 
tice that the difference [f$ eq — fQ eq ]y = _u 2 can De evaluated in two ways. On one hand, 
[fr\=-L/2 = [fi] y =-L/2 ~ [ft\=-L/2, where the first contribution is due to ([22]), and 
the second is due to the stationary solution for the equilibrium ([2]) with the velocity ([TBI) , 
whereupon [/ 5 ncq - fT\ = _ L/2 = ^(Ui+^(-^+rV)). On other hand, using ([21]), we find 
[f^ ~ /6 eq ]j/=-L/2 = (if - Matching these expressions, we find a relation between integration 
constants n and V: 



71 



(23) 



3+^5 (i-rvy 

Similarly, at the top plate (y = L/2, see Fig. [[]), / 4)7 8 |^ = _ L/2 = f^ 7>a (p, pU 2 , 0). Again, 
computing the difference [/° eq — f% e(V \ y=L /2 i n ^ w0 ways as described above, we find: 

V3c s U 2 



71 



3 + 7? {% + tV) 



(24) 



Comparing (1231) and (1241) . we find coefficients V and it, and, making use of ([18]) . we arrive 
at the result for the velocity profile ([!]), while the non-equilibrium shear stress P£y q = pvr 
reads 

where Og is given by ([5]). Note that in the D2Q1Q model the result for the shear stress 



is obtained by replacing 9 with 16 ([7]). Results for the shear stress are compared with 
the data of Willis [l5j and DSMC simulation in Fig. [31 In Tab. [H the limiting values of 
the effective shear viscosity, v e g = — P^ q L(A[/p) _1 , at the infinite Knudsen number for the 
D2Q9 and the D2Q16 LB models are compared with the Boltzmann-BGK result flo) ]. 
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FIG. 3: Left: Shear stress at various Knudsen numbers. Labels are as in Fig. [2j Plotted is the 
reduced function P* y = P xy /P^ where P££ is the shear stress at Kn — > oo of the Boltzmann-BGK 
model lj]]. Right: Nonequilibrium normal stress difference at Kn = 0.6. Line: solution (I26j) : 
Symbol: DSMC simulation. 



D2Q9 


0.723 


D2Q16 


1.113 


BGK [15] 
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TABLE I: Effective shear viscosity v g r at Kn — ► oo, in units of 



Same method is used in order to evaluate the two remaining integration constants A and 
B. Namely, we evaluate /^ eq (bottom plate) and f£ eq (top plate) in the two ways as de- 
scribed above. After some algebra, we find the results for the non-equilibrium normal stress 
difference and the transversal energy flux 

/ AU\ 2 



ncq 



TV 



P 



V L J (1 

2 



-P 



AU 



2Kn) 2 
v 



2 — e 2Kn cosh 



y 

KnLJl ' 



(26) 



2Kn) s 



2y — KnLe 2Kn sinh 



y 



KnL 



+2UF%*. (27) 



Expressions for the velocity u x (jlj), the shear stress P^ q ( 1251) . the normal stress difference 
A^ neq (1261) , and the energy flux ([271), when substituted into fa = /f q + /f eq (see (J2D and 
(pTj) ). furnish the exact solution of Couette flow for the nonlinear D2Q9 model. The same 
solution method is applicable to any member of the LB hierarchy although algebra becomes 
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more involved. Solution of the D2Q16 model leading to the result ([6]) will be presented in 
a detailed publication. 

We have already mentioned that the solution for the velocity (jl]) is Galilean invariant. 
Same holds also for the higher-order moments. Indeed, the stress tensor is obviously Galilean 
invariant (that is, -P°^ q (1251) and iV neq (125]) depend only on [7 2 — U±), and also the energy 
flux transforms correctly under U[ 2 — ^1,2 + C: from the definition of the energy flux 
l' y — % + 2CP xy + C 2 j y , which is manifested in ( 127]) (j y = in the present solution). The 
normal stress difference (|26|) is a positive-definite function which is consistent with kinetic 
theory of gases. Importantly, the nonvanishing of iV neq and g° eq is the direct implication of 
the nonlinearity of the kinetic equation ([T]) (manifested by the AU 2 dependence in (I2b'p and 
(12"?]) ). and cannot be predicted on the basis of linearized kinetic theory Q, [3]. For that 
reason, the DSMC method [3] was used in order to validate the solution in the nonlinear 
domain. Parameters of the DSMC simulation correspond to the hard sphere model of Argon 
at normal conditions and diffuse scattering at the walls was implemented. The value of 
the velocity U\ = —U2 = 0.5C S , where C s = ^J^k^T /3m is the speed of sound of the 
three-dimensional ideal gas, was used to maintain a quasi-isothermal flow. The normal 
stress difference (1261) is mapped onto DSMC data in Fig. [3] which qualitatively confirms the 
prediction. 

Finally, in view of the concerns mentioned in the introduction, we validate the lattice 
Boltzmann method for the kinetic equation ([I]). The LB space-time discretization of ([1]) is 
standard, see, e. et. Simulation setup used corresponds to U\ = 0, periodic boundary 
conditions were applied in the x-direction. Simulation were run till the steady state was 
reached, a grid convergence study was also performed. It is evident in Fig. H] that the 
simulated relative slip accurately reproduces the exact solution (jl]). Same holds for all other 
moments. Thus, applications of the LBM to microflows are by no means an "artifact of 



numerics" [12j]. Importance of physically relevant boundary conditions must be stressed. 
We have verified that if the standard "bounce-back" boundary condition of the LB method 
is used in the present problem, then the analytical result predicts vanishing slip velocity 
at all Knudsen numbers; thus, results of micro-flow simulations with the "bounce-back" 
boundary conditions or their derivatives are questionable indeed [la ]. 
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FIG. 4: Exact solution for the relative slip velocity (line) and the lattice Boltzmann simulation 
(symbol) function of Knudsen number. 

To conclude, our analytical results suggest that the hierarchy of lattice Boltzmann models 
is the way to approximate the kinetic theory. Without a denial of a body of LB simulations, 
it must be appreciated that only the exact solutions answer unambiguously the question of 
the physical validity of the method. Our results reveal that applications to LB methods to 
microflows should be based on LB models with larger velocity sets if one seeks a quantitative 
prediction, especially in the transient regime. It would be quite interesting to extend the ap- 
proach developed herein to obtain closed- form solutions in other cases of kinetic theory [14j , 
and eventually also to three-dimensional cases. In that respect, we were able to extend the 
present analysis to the three-dimensional D3Q27 model, results will be reported elsewhere. 

We acknowledge useful discussions of some of the results with G. Doolen, H.C Ottinger, 
S. Succi and V. Yakhot. Support by BFE Project 100862 (I.V.K.), by the ETH Project 
0-20235-05 (N.I.P), by CCEM-CH (I.V.K. and S.A.) and by NTU Office of Research (S.A.) 
is gratefully acknowledged. 
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